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Abstract 

Background: Functional studies liave demonstrated tliat microRNAs (miRNAs or miRs) play critical roles in a wide 
spectrum of biological processes including development and disease pathogenesis. To investigate the functional 
roles that miRNAs play during chicken skeletal muscle development, the miRNA transcriptomes of skeletal muscles 
from broiler and layer chickens were profiled using Solexa deep sequencing. 

Results: Some miRNAs have multiple isoforms and several miRNAs^ are present at higher levels than their 
corresponding miRNAs. Thirty three novel and 189 known chicken miRNAs were identified using computational 
approaches. Subsequent miRNA transcriptome comparisons and real-time PGR validation experiments revealed 17 
miRNAs that were differentially expressed between broilers and layers, and a number of targets of these miRNAs 
have been implicated in myogenesis regulation. Using integrative miRNA target-prediction and network-analysis 
approaches an interaction network of differentially expressed and muscle-related miRNAs and their putative targets 
was constructed, and miRNAs that could contribute to the divergent muscle growth of broiler and layer chickens 
by targeting the ACVR2B gene were identified, which can causes dramatic increases in muscle mass. 

Conclusions: The present study provides the first transcriptome profiling-based evaluation of miRNA function 
during skeletal muscle development in chicken. Systematic predictions aided the identification of potential miRNAs 
and their targets, which could contribute to divergent muscle growth in broiler and layer chickens. Furthermore, 
these predictions generated information that can be utilized in further research investigating the involvement of 
interaction networks, containing miRNAs and their targets, in the regulation of muscle development. 



Background 

Embryonic patterning and organogenesis involve coordi- 
nated differentiation, migration, proliferation and pro- 
grammed cell death in metazoans. These complex 
cellular and developmental processes rely on precise 
spatiotemporal networks that regulate transcription fac- 
tors at multiple levels including mRNA transcription 
and translation, protein stability and degradation. 
Recently, evidence has demonstrated that microRNAs 
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(miRNAs or miRs) are involved in diverse aspects of 
biology including developmental regulation and the 
pathogenesis of human diseases [1-4]. miRNAs are small 
19-24 nucleotide (nt) regulatory RNAs that generally 
modulate gene expression through translational repres- 
sion or by causing deadenylation and degradation of tar- 
get mRNAs [5,6]. However, miRNAs could function as 
activators to regulate gene expression [7,8]. The biogen- 
esis of miRNAs is spatiotemporally regulated by various 
mechanisms [9], providing additional evidence that miR- 
NAs are functionally significant, and potentially key reg- 
ulators of gene expression during development [10-15]. 

An essential role for miRNAs in terms of regulating 
skeletal muscle development is evident from studies 
demonstrating that deletion of a conditional Dicer allele 
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in embryonic skeletal muscle results in perinatal lethal- 
ity due to skeletal muscle hypoplasia [16]. In particular, 
the critical roles of three muscle-specific miRNAs, miR- 
1, miR-133 and miR-206, in the regulation of myogen- 
esis have been well documented [17,18]. miR-1 and 
miR-133 have been reported to regulate different aspects 
of skeletal muscle development in vitro and in vivo [19]. 
miR-1 promotes myocyte differentiation by repressing 
the expression of histone deacetylase 4 (HDAC4), a 
negative regulator of differentiation and a repressor of 
the MEF2 (myocyte enhancer factor-2) transcription fac- 
tor [19]. In C2C12 myoblasts, miR-133a promotes prolif- 
eration, in part, by repressing serum response factor 
(SRF) [19]. Like miR-1, miR-206 promotes differentia- 
tion of C2C12 myoblasts in vitro, miR-206 induces mus- 
cle differentiation by repressing the expression of the 
DNA polymerase a subunit (Polal) [20], connexin 43 
(Cx43) [21], follistatin-like 1 (FstU) and utrophin (Utrn) 
[22]. In addition to muscle-specific miRNAs, several ubi- 
quitously expressed miRNAs have a role to play during 
muscle development. For example, zebrafish miR-214 
was reported to regulate the slow muscle phenotype by 
targeting suppressor of fused {Sufu), a negative regulator 
of hedgehog signaling [23]. Expression of miR-181 iso- 
forms, miR-lSla and miR-181b, are induced upon initia- 
tion of myogenesis and they participate in the regulation 
of myoblast differentiation by repressing HoxA-11 pro- 
tein levels [24]. The functional significance of miRNAs 
in terms of controlling myogenesis has been documen- 
ted, but the majority of miRNAs are abundantly 
expressed. Therefore, identif)^ing novel miRNAs that are 
expressed at low levels during skeletal muscle develop- 
ment but are functionally important requires robust 
approaches such as high-throughput deep sequencing 
technology. 

The chicken {Gallus gallus) is an established model 
organism for studying vertebrate development, primarily 
because chicken embryos are readily accessible and 
easily manipulated [25]. In addition, a variety of stan- 
dard chicken breeds with different phenotypes are read- 
ily available, which collectively represent a valuable 
genetic resource. Broiler chickens (bred for meat pro- 
duction) and layer chickens (bred for egg production) 
are ideal model systems for studying the molecular 
mechanisms underlying myogenesis [26]. During the 
past 80 years, genetic selection in broilers has concen- 
trated on a high growth rate and large muscle mass; in 
contrast, layers have been selected for egg production. 
Therefore, even under optimal growth conditions, the 
body size of layers is smaller than that of broilers owing 
to intrinsic genetic differences between the two varieties. 
These unique biological features of broilers and layers 
allow muscle development to be investigated. In pre- 
vious studies, we have successfully identified protein- 



coding and non-coding genes with roles during myogen- 
esis using broilers and layers as model systems [27,28]. 

In the current study, previous work was expanded to 
identify miRNAs involved in myogenesis regulation by 
comparing the miRNAs transcriptome in skeletal muscle 
tissues of broilers and layers. Solexa deep sequencing was 
carried out to profile miRNAs expressed in chicken ske- 
letal muscle tissues. Sequence-tag analyses have shown 
that a group of highly abundant, known miRNAs are 
expressed in skeletal muscles and 33 novel putative 
chicken miRNAs from skeletal muscle tissues have been 
identified. Comparing the expression patterns of known 
and novel miRNAs demonstrated that they were signifi- 
cantly differentially expressed between broiler and layer 
chicken muscle tissues. These results were confirmed 
using microarrays and real-time reverse transcription- 
polymerase chain reaction (RT-PCR) validation experi- 
ments. Of the 17 miRNAs examined using RT-PCR, nine 
presented with an expression pattern consistent with the 
microarray analysis; 15 miRNAs had a pattern consistent 
with the deep sequencing data. Using computational pre- 
diction, targets for these differentially expressed miRNAs 
and muscle-related miRNAs were identified, and an 
interaction network was constructed. Furthermore, miR- 
1 was demonstrated specifically to target the 3' untrans- 
lated region of the activin receptor IIB gene, ACVR2B, 
which can cause dramatic increases in muscle mass [29]. 
This integrative analysis highlights the complexity of 
gene expression networks regulated by microRNAs in 
muscle cells during muscle development. 

Results 

Characterization of the miRNA transcriptome of skeletal 
muscle from broiler and layer chickens using deep 
sequencing 

Solexa sequencing was used to profile miRNAs 
expressed in layer and broiler chicken skeletal muscles. 
Sequencing of a small RNA fraction (16-30 nt) from 
total RNA extracted from pectoralis muscles collected 
from 10-day-old chicken embryos yielded 2,700,003 and 
2,576,562 reads for the layer and broiler libraries, 
respectively (Figure lA). Of these, 1,987,912 layer 
sequences and 1,553,308 broiler sequences, which 
account for more than 67% of the total reads, were per- 
fectly mapped to the chicken genome (May 2006). The 
sequencing data were simplified by grouping all identical 
sequence reads together; therefore, 105,475 unique layer 
sequences and 89,148 unique broiler sequences were 
used for subsequent analysis (Figure lA). 

The most abundant size class in the small RNA 
sequences distribution was 22 nt, followed by 21 and 23 
nt (Figure IB and IC), and this was consistent with the 
known 21-23 nt range for miRNAs. To assess the effi- 
ciency of deep sequencing for miRNA detection, all 
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Figure 1 Deep sequencing results and annotations of small RNAs from chicken skeletal muscle. A. Number of small RNA reads from 
broilers and layers. 8, C Size distribution of sequenced small RNAs. D. Annotations of sequenced small RNAs. 




sequence reads were annotated and classified by analyz- 
ing the sequence tags in relation to the data from miR- 
Base (version 16), RefSeq mRNA, RepeatMasker and 
non-protein-coding RNAs annotated by ENSEMBL. The 
sequence tag annotations demonstrated that known 
chicken miRNAs (gga_miRNAs) and metazoan miRNA 



homologs accounted for -50% of all sequence reads in 
the broiler and layer libraries (Figure ID). These results 
indicate that the deep sequencing data were highly 
enriched for mature miRNA sequences, suggesting that 
the data are reliable for expression profiling of known 
miRNAs and deep mining for novel miRNAs. 
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To investigate the expression of known miRNAs in 
broiler and layer skeletal muscles, the numbers and dis- 
tribution of small RNA sequences that matched known 
chicken miRNA genes were analyzed. The results 
demonstrated that of 467 known chicken miRNAs and 
77 miRNA^'s in the miRBase (version 16), perfect 
matches to 231 miRNAs and 29 miRNA^'s were obtained 
in the sequencing data (Figure 2 A and Additional file 1). 
Among the sequences that were not perfectly matched 
to known chicken miRNAs or miRNA ''s there were 244 
metazoan miRNA homologs and 72 metazoan miRNA'' 
homologs (Figure 2 A and Additional file 2). 
As presented in Additional file 5 and 6, known miRNAs 
and metazoan homologs had a broad range of 



expression levels in skeletal muscle tissues, ranging from 
hundreds of thousands of sequence reads for the most 
abundant miRNAs to single reads for the least abun- 
dant The distribution of read numbers for the known 
miRNAs is summarized in Figure 2B, The 33 most 
abundant miRNAs (i.e. those with > 1,000 reads) are 
presented in Table 1. 

miRNA transcriptome analysis demonstrated the presence 
of several highly abundant miRNAs in skeletal muscles of 
broilers and layers 

Almost all known muscle-specific miRNAs (myomiRs) 
were represented among the more abundant miRNAs 
identified. The most abundant miRNA was gga-miR- 
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Figure 2 Known miRNAs and homologs of metazoan miRNAs detected in chicken skeletal muscle. A. Numbers of known miRNAs and 
chicken liomologs of metazoan miRNAs detected by perfect matcli in tine present study. B. Reads distribution of known miRNAs and cliicken 
liomologs of metazoan miRNAs in cliicken skeletal muscle. 



Li et al. BMC Genomics 201 1, 12:186 
http://www.biomedcentral.eom/1 471-21 64/1 2/1 86 



Page 5 of 20 



Table 1 The most abundant miRNAs in chicken skeletal 
muscles as determined using deep sequencing 



miRNA ID 


Lay6r 


Broil6r 
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^tudip^ rpl;)tprl tn ^kplptril 
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131609 


0.59 


[20-22] 


gga-let-7c 


32663 


87111 


2.67 
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206, which was represented by approximately 200,000 
sequence reads in the broiler and layer libraries (Table 
1). The predominance of miR-206 is consistent with its 
well established function during skeletal muscle develop- 
ment [30] and reported role during chicken myogenesis 
[31,32]. Two other myomiRs, miR-1 [19] and miR-181 
[24], were high-count sequences in both libraries (Table 
1). Compared with the three myomiRs, forms of the 
myomiR miR-133 were expressed at low levels in the 
skeletal muscle libraries: there were 126 reads for gga- 
miR-133a in the layers library and 67 in the broilers 
library; 7 reads for gga-miR-133b in the layers library 
and 12 in the broilers library; one gga-miR-133c read in 
the layers library and no read in the broilers library 



(Additional file 1). These variations in abundance could 
reflect differences in the roles of these miRs in terms of 
the regulation of myogenesis [19,24]. In addition to 
miR-206, miR-1 and miR-181, nine other miRNAs 
among the most abundant in these libraries (miR-221, 
miR-222, miR-21, miR-103, miR-130, miR-99, miR-30, 
miR20, and miR128) have been implicated in the prolif- 
eration and differentiation of muscle cells (Table 1) 
[15,19,33]. Therefore, the miRNA transcriptome for 
layers and broilers revealed by this analysis is highly 
enriched for miRNAs involved in myogenesis regulation. 

The expression patterns of several mlRNA^s during 
chicken siceletal muscle development are unique 

A total of 29 miRNA ''s were detected in broiler and 
layer libraries (Figure 2A and 2B). The majority were 
expressed at low levels (Additional file 1), but of those 
expressed at higher levels, such as miR-lSla and miR- 
1677, the read counts were significantly lower than 
those of the corresponding miRNAs (Table 2). One 
striking exception to this general trend was gga-miR- 
140''; this was present as 911 reads in the layers and 711 
reads in the broilers libraries, but gga-miR-140 was not 
detected in either library (Table 2), suggesting that gga- 
miR-140'' functions during chicken skeletal muscle 
development. gga-miR-126'' is another case in which 
only miRNA'' was detected (Table 2). There were several 
cases, such as miR-199 and miR-1329, where miRNA 
and miRNA" were generated at similar levels (Table 2). 

Analysis of sequence variants indicated that many 
miRNAs possess isomiRs 

As found in previous deep sequencing studies, heteroge- 
neity at the 5' and/or 3' ends of miRNAs was observed 
(Figure 3, Additional file 3). miRNAs with such varia- 
tions from their miRBase reference sequences are 
referred to as isomiRs [34,35]; some typical examples 
are presented in Figure 3. In the majority of cases (e.g. 
gga-miR-221) the most abundant isoform is identical to 
the reference in miRBase (Figure 3). In some cases, such 
as gga-miR-222 and gga-miR-128, more than one highly 
abundant isoform was present (Figure 3), indicating that 

Table 2 A comparison of read counts between miRNA 
and the corresponding miRNA^ 



miRNA ID miRNA miRNA^ 





Layer 


Broiler 
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gga-miR-1 81 a 
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2543 


142 
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gga-miR-1 329 


16 


123 


36 


44 


gga-miR-140 


0 


0 
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Broiler 
0 



>eea-miR-221 

ACTGTCCAGGTTTGGGGCATGAACCTGGCATACAATGTAGATTTCTGTGTTTGTTAAGCAACAGCTACATTGTCTGCTGGGTTTCCAGCTGCCTGGAA 

....(((((((....((...(((((..(((.(((((((((....((((((( ))).))))))))))))).)))..)))))))....))))))). 

CTACATTGTCTGCTGGGTTTC 

GCTACATTGTCTGCTGGGTT 

GCTACATTGTCTGCTGGGTTT 

GCTACATTGTCTGCTGGGTTTC 

AGCTACATTGTCTGCTGGG 

AGCTACATTGTCTGCTGGGT 

AGCTACATTGTCTGCTGGGTT 

AGCTACATTGTCTGCTGGGTTT 

AGCTACATTGTCTGCTGGGTTTC * 

AGCTACATTGTCTGCTGGGTTTCC 

AGCTACATTGTCTGCTGGGTTTCCA 

CAGCTACATTGTCTGCTGGGTTTC 

>eea-miR-222 

GTAGTTGCCCATCAATCGCTCAGTAGTCAGTGTAGATTCTGTCTTTACAATCAGCAGCTACATCTGGCTACTGGGTCTCTGATGACATCTCATATCT 

(.((.((..(((((...((((((((((((((((((...((( )))...))))).)))))))))))))...))))).)).))) 

CTACATCTGGCTACTGGGTCTC 

GCTACATCTGGCTACTGGGTCT 

GCTACATCTGGCTACTGGGTCTC 

GCTACATCTGGCTACTGGGTCTCT 

AGCTACATCTGGCTACTGG 

AGCTACATCTGGCTACTGGG 

AGCTACATCTGGCTACTGGGT 

AGCTACATCTGGCTACTGGGTC 

AGCTACATCTGGCTACTGGGTCT 

AGCTACATCTGGCTACTGGGTCTC * 

AGCTACATCTGGCTACTGGGTCTCT 

CAGCTACATCTGGCTACTGGGTCT 

CAGCTACATCTGGCTACTGGGTCTCT 

>eea-miR-128-l 

GAGCTGTTGAATTCGGGGCCGTAACACTGTCTGAGAGGTTTATATTTCTCACAGTGAACCGGTCTCTTTTTCAGCTGCTTC 
((((.((((((...(((((((...((((((..(((((((....)))))))))))))...)))))))...)))))).)))). 

ACAGTGAACCGGTCTCTT 

ACAGTGAACCGGTCTCTTT 

ACAGTGAACCGGTCTCTTTT 

ACAGTGAACCGGTCTCTTTTT 

CACAGTGAACCGGTCTCTT 

CACAGTGAACCGGTCTCTTT 

CACAGTGAACCGGTCTCTTTT 

CACAGTGAACCGGTCTCTTTTT 

CTCACAGTGAACCGGTCTCTTT 

TCACAGTGAACCGGTCTCT 

TCACAGTGAACCGGTCTCTT 

TCACAGTGAACCGGTCTCTTT * 

TCACAGTGAACCGGTCTCTTTT 

TCACAGTGAACCGGTCTCTTTTT 

TCACAGTGAACCGGTCTCTTTTTC 

TCTCACAGTGAACCGGTCTCT 

TCTCACAGTGAACCGGTCTCTTT 

>2ea-miR-181a-l 

GTAGTGGTTGCTTCAGTGAACATTCAACGCTGTCGGTGAGTTTGGAATTTAAGTGAAAACCATCGACCGTTGATTGTACCCTCCAGCTAACCATCCTCCTCCT 

...(((((((((..((.(.(((.((((((..(((((((.((((...((....))..))))))))))))))))).))).).))..))).)))))) 

CATTCAACGCTGTCGGTGAGTTT 

ATTCAACGCTGTCGGTGAGTT 

ATTCAACGCTGTCGGTGAGTTT 

ACATTCAACGCTGTCGGTGAG 

ACATTCAACGCTGTCGGTGAGT 

ACATTCAACGCTGTCGGTGAGTT 

ACATTCAACGCTGTCGGTGAGTTT 

AACATTCAACGCTGTCGGTGA 

AACATTCAACGCTGTCGGTGAG 

AACATTCAACGCTGTCGGTGAGT * 

AACATTCAACGCTGTCGGTGAGTT 

AACATTCAACGCTGTCGGTGAGTTT 

AACATTCAACGCTGTCGGTGAGTTTG 

GAACATTCAACGCTGTCGGTGA 

G AACATTCAACGCTGTCGGTGAG 

GAACATTCAACGCTGTCGGTGAGT 

GAACATTCAACGCTGTCGGTGAGTT 

GAACATTCAACGCTGTCGGTGAGTTT 

>eea-miR-la-l 

GAGAGACATACTTCTTTATATGCCCATATGAACCTGGCAATCTATGGAATGTAAAGAAGTATGTATTTCA 

((((.((((((((((((((((..(((((.(( ))))))).)))))))))))))))).)))). 

GAATGTAAAGAAGTATGT. . 

GAATGTAAAGAAGTATGTAT. . . 

GAATGTAAAGAAGTATGTATT. . 

GAATGTAAAGAAGTATGTATTT. . 

GGAATGTAAAGAAGTATGT. . . 

GGAATGTAAAGAAGTATGTA. . 

GGAATGTAAAGAAGTATGTAT. . . 

GGAATGTAAAGAAGTATGTATT. . 

TGGAATGTAAAGAAGTATG. . . 

TGGAATGTAAAGAAGTATGT. , 

TGGAATGTAAAGAAGTATGTA. . . . 

TGGAATGTAAAGAAGTATGTAT. . . 

TGGAATGTAAAGAAGTATGTATT. . 

TGGAATGTAAAGAAGTATGTATTT. . 

ATGGAATGTAAAGAAGTATGT. . . 

ATGGAATGTAAAGAAGTATGTA. . 

ATGGAATGTAAAGAAGTATGTAT. . . 

TATGGAATGTAAAGAAGTATGT. . 

>2ea-miR-499 

TTTGAGGGAGCGGCAGTTAAGACTTGTAGTGATGTTTAGATAATGTATTACATGAACATCACTTTAAGTCTGTGCTACTTCTCTCCTCAT 
..((((((((.((.(((..(((((((.(((((((((((..(((....)))..))))))))))).)))))))..))).)).))).))))). 

AAGACTTGTAGTGATGTTT 

AAGACTTGTAGTGATGTTTA 

AAGACTTGTAGTGATGTTTAGA 

TAAGACTTGTAGTGATGTTT 

TAAGACTTGTAGTGATGTTTA 

TAAGACTTGTAGTGATGTTTAGA 

TTAAGACTTGTAGTGATGTTT 

TTAAGACTTGTAGTGATGTTTA 

TTAAGACTTGTAGTGATGTTTAG * 

Figure 3 IsomiRs from several gga-miRs. Reads alignments of the various isoforms of several gga-miRs are presented. The sequence of the 
gga-miR hairpin is presented in the top line; the brackets below denote the secondary structure. Reads that aligned with the mature gga-miR 
sequence as reported in miRBase are denoted by a series of asterisks. The number of reads corresponding to each sequence is presented on the 
right. 
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some miRNAs have more than one functional isoform 
in specific tissues/organs. For some miRNAs, such as 
gga-miR-181a, gga-miR-la and gga-miR-499, the most 
abundant isoform was not among the known miRNA 
sequences reported in miRBase 16 (Figure 3). A similar 
phenomenon has previously been observed for miRNAs 
identified in chicken embryos [36], suggesting that a 
refinement to the miRBase annotations for chicken miR- 
NAs is required to reflect experimentally observed 
abundances. 

Novel miRNAs are less abundant and less evolutionarily 
conserved in chicken skeletal muscle 

In addition to profiling known miRNAs, deep sequen- 
cing is a powerful strategy for discovering novel miR- 
NAs that may not have been detected using traditional 



methods for sequencing cDNA libraries. Using the miR- 
Deep program as a predictive tool [37], 33 putative 
novel chicken miRNAs were obtained from broiler and 
layer sequence tags (Table 3, Additional file 4). Genomic 
sequence analyses demonstrated that six of these puta- 
tive miRNAs were located in the exons of annotated 
genes, fourteen resided in the introns of annotated 
genes and thirteen were present in intergenic regions 
(Additional file 5). The putative novel miRNAs were less 
abundant (Table 3) than known miRNAs (Figure 2A 
and 2B). Only one novel putative miRNA had read 
counts greater than 100 in the library (Table 3). 

To investigate evolutionary conservation of the 33 
novel chicken miRNAs, a search for highly similar 
sequences among human, mouse, rat, opossum, frog and 
zebrafish genomic sequences was carried out using a 



Table 3 Novel chicken miRNAs predicted by mIRDeep 



miRNA ID 


chr 


strand 


start 


end 


Layer^ 


Broiler^ 


gga-miR-Nl 


chr3 


+ 


434S337 


434S3SS 


13S 


62 


gga-nniR-N2 


clirUn_random 


- 


13277374 


1327739S 


2 


S2 


gga-miR-NS 


clirUn_random 


+ 


30704197 


3070421 S 


29 


21 


gga-miR-N4 


chrl7 


+ 


247706S 


24770S7 


18 


14 


gga-miR-NS 


chr27 


+ 


1117232 


11172SS 


22 


6 


gga-nniR-N6 


chr2 


+ 


1030429S7 


103042979 


26 


0 


gga-miR-N7 


Chris 


- 


S042929 


S0429S2 


20 


S 


gga-miR-NS 


clir4 




3214463 


32144S4 


23 


0 


gga-miR-N9 


chr6 


+ 


31SS03S7 


31SS037S 


18 


4 


gga-miR-Nl 0 


chrl 




1S021S042 


1S021S06S 


14 


3 


gga-miR-Nl 1 


clir4 




S12S7SS1 


S12S7S72 


7 


9 


gga-miR-Nl 2 


chrZ 




2S24721S 


2S24723S 


12 


4 


gga-miR-Nl 3 


chr3 


+ 


SSI 17942 


SS11796S 


13 


1 


gga-miR-N14 


clir4 


+ 


S633S2 


S63373 


11 


3 


gga-miR-Nl 5 


chr26 




2669S12 


2669S34 


11 


2 


gga-miR-Nl 6 


chr3 


+ 


S930796S 


S9307990 


5 


7 


gga-miR-Nl 7 


Chris 




41162SS 


4116279 


6 


S 


gga-miR-N18 


chr3 


+ 


724977 


724994 


0 


7 


gga-miR-N19 


chr4 


+ 


21S123S 


21S1261 


4 


3 


gga-miR-N20 


chr7 




37S93S31 


37S93SS1 


6 


1 


gga-miR-N21 


chrl 9 


+ 


4S62173 


4S62194 


1 


4 


gga-miR-N22 


chr20 




10S96392 


10S96413 


5 


0 


gga-miR-N23 


Chris 




112SS901 


112SS922 


4 


0 


gga-miR-N24 


chr24 


+ 


261 SSI S 


261SS40 


3 


1 


gga-miR-N25 


chr27 




4426392 


4426414 


4 


0 


gga-miR-N26 


chr4 


+ 


1662S71S 


1662S736 


3 




gga-miR-N27 


chrUn_random 




402S0702 


402S072S 


3 




gga-miR-N28 


chrl 


+ 


S2701699 


S2701723 


2 




gga-miR-N29 


chrlO 




1641S1S1 


1641S173 


1 




gga-miR-N30 


chr2 


+ 


1333037S6 


133303776 


1 




gga-miR-N31 


chr27 




39S7S24 


39S7S4S 


2 


0 


gga-miR-N32 


chr3 




49S36173 


49S36194 


1 


1 


gga-miR-N33 


chr7 


+ 


12SS0964 


12SS09SS 


0 


2 



*Sequenced read-numbers are presented. 
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BLAST analysis. Obtaining mature miRNAs from 
homology sequences does not necessarily signify that 
the miRNAs are conserved as they might not be capable 
of forming hairpin structures. We further identified 
hairpin-like RNA structures using RNAfold (see Materi- 
als and Methods). The same analysis was carried out for 
known chicken miRNAs. The results demonstrated that 
novel miRNAs are less evolutionarily conserved (Addi- 
tional file 6), a result that is consistent with previous 
studies [36]. Further analyses using a multiple alignment 
of six vertebrate genomes in the UCSC database with 
that of the chicken genome revealed that only gga-miR- 
N3 was conserved in at least one of the analyzed verte- 
brate genomes. gga-miR-N3 exists in zebrafish and frogs 
but is lost after the emergence of mammalian lineages 
(Additional file 7). The remaining 32 miRNAs could be 
avian- and/or chicken-specific miRNAs. To identify 
potential chicken-specific miRNAs, the sequences were 
checked against the Zebra Finch {Taeniopygia guttata) 
Alignment Net Track in UCSC. Nineteen of the 32 
novel chicken miRNAs were present in the Zebra Finch, 
suggesting that the remaining 13 novel miRNAs could 
be specific to the chicken lineage (Additional file 8). 
Combining conservation and relative abundance infor- 
mation for the newly identified and known miRNAs 
revealed that the evolutionarily conserved miRNAs were 
among the most abundant, supporting a correlation 
between evolutionary conservation and the expression 
level of miRNAs. 

Identification of differentially expressed miRNAs in broiler 
and layer skeletal muscle 

The main objective of the present study was to identify 
miRNAs involved during skeletal muscle development 
by comparing skeletal muscle miRNA transcriptomes in 
broilers and layers. Analysis of sequencing results 
demonstrated that more than 80% of reads overlapped 
between broilers and layers (Figure 4A). The overlap 
between libraries was greater (94%) for those reads with 
perfect genomic matches (Figure 4B), suggesting that 
the deep sequencing data were reliable for direct com- 
parison of miRNA abundance between broilers and 
layers. 

In addition to the 33 novel miRNAs, 189 known miR- 
NAs were identified in miRBase using miRDeep (Addi- 
tional file 9). Comparing Table S6 with Table SI 
demonstrates that the number of reads could differ for 
the same miRNA. In Table SI, the reads number for 
each miRNA is based on perfect matches to known 
chicken miRNAs in miRBase. As presented in Figure 3 
and Additional file 3, many miRNAs have different iso- 
miRs in addition to perfect matches. Counting only per- 
fect-match isoforms may not be appropriate as the 
isoform listed in the miRBase may not be the only 



functional isoform. In miRDeep, different isoforms of 
the same miRNA are counted together. To arrive at the 
figure of 189 known miRNAs identified by miRDeep, 
the DEGseq package [38] was used to identify differen- 
tially expressed miRNAs on the basis of potentially sig- 
nificant changes in relative miRNA abundance between 
broilers and layers. Expression of 102 miRNAs was sig- 
nificantly different between broilers and layers (Addi- 
tional file 10, Figure 4C). miRNA microarrays were 
employed to characterize the expression profiles of these 
102 differentially expressed miRNAs further (Additional 
file 10, Figure 4D). 

To validate the differential expression of these miR- 
NAs between broilers and layers, 17 miRNAs were ran- 
domly selected and their expression levels quantified 
using real-time RT-PCR (Figure 5). Of the 17 miRNAs 
examined, nine (52.9%) had an expression pattern con- 
sistent with the microarray analysis (Figure 5A and 5C), 
and 15 miRNAs (88.2%) presented with a pattern con- 
sistent with the deep sequencing data (Figure 5A and 
5B). These data provide evidence that deep sequencing 
is a more sensitive and reliable method for identifying 
differentially expressed miRNAs than miRNA 
microarrays. 

Target prediction and network analysis highlight the 
complexity of interactions among miRNAs and their 
targets during muscle development 

Of the 17 differentially expressed miRNAs confirmed 
using real-time RT-PCR, six have been functionally 
linked to myogenesis [39-42]. However, the majority 
including one novel miRNAs (gga-miR-N2) and eight 
known miRNAs (miR-101, miR-15b, miR-15c, miR- 
1677, miR-200, miR-460, gga-mir-2188 and miR-429) 
have not been implicated in the regulation of muscle 
development. To approach the question of how miRNAs 
could function in concert with their target genes in 
terms of controlling muscle development and to provide 
some molecular insight into the process, targets of the 
miRNAs were identified and a possible regulatory net- 
work of interactions among miRNAs and their targets 
was constructed. The strategy and workflow are sum- 
marized in Figure 6A. 

The starting point of the miRNA target prediction 
strategy was the 16 validated, differentially expressed, 
known miRNAs and eight other muscle-related miR- 
NAs. TargetScan (version 5.1) [43] was used to predict 
the putative targets for these 24 miRNAs, identifying 
more than 1000 annotated mRNA transcripts that were 
potential targets (Additional file 11). The mechanism of 
miRNA function predicts that miRNAs and their targets 
normally exhibit correlated expression patterns [44]. 
Therefore, an mRNA transcriptome analysis was per- 
formed using microarrays to identify mRNAs in 



Li et al. BMC Genomics 201 1, 12:186 Page 9 of 20 

http://www.biomedcentral.eom/1 471 -21 64/1 2/1 86 



A 



D 



Layer Broiler 



h4 I-? 1-4 CQ CQ CQ 




348129 4241870 1686566 



B 



Layer Broiler 




gga 

gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 
gga 



-iniR-101 

-miR-30a 
-miR-456 
-miR-30d 
-miR-1456 
-miR-31 
-rtiiR-la-1 
-miR-200a 
-miR-206 
-miR-lOa 
-miR-383 
-miR-101-2 
-miR-133a-l 
■miR-181b-2 
let-7i 
-miR-la-2 
■miR-181a-l 
-miR-203 
-miR-33-1 
-miR-1557 
-miR- 103-2 
-miR-181a-2 
miR-lSlb-l 
-niiR-30e 
miR- 107 
-let-7f 
miR- 167 7 
-miR- 193a 
-miR-184 
miR- 168 4 
-miR-365-1 
-miR-30c-l 
-miR-103-1 
miR- 172 8 
-mir-3532 
-miR- 14 8a 
-miR- 1671 
-miR-133a-2 
-miR-1716 
-miR- 14 6c 
-miR- 10b 
-miR- 15 60 
-miR-138-2 
-miR-138-1 
miR- 16 98 
-miR-34b 
■let-7a-2 
-miR-30b 
-miR- 4 60 
-miR- 18b 
-let-7a-l 
-miR-130c 
-miR- 16-1 
-miR- 2 3b 
-miR-1416 
-miR-128-1 
-miR-106 
-miR-196-3 
■miR- 454 
-miR- 15b 
-let-7c 
-miR-130a 
-miR-24 
-let-7a-3 
-miR- 13 97 
-miR- 18a 
-miR-15G 
-let-7j 
-miR- 100 
-miR-222 
-miR-21 
miR-196-1 
-mir-2188 
-miR-199-1 
-miR- 48 9 
miR- 9 9a 
miR- 2 00b 
-miR-205b 
-miR- 16 -2 
-miR-21 5 
-miR-7-3 
-miR-551 
-miR- 42 9 
let-7g 
-miR- 92 
-miR- 2 7b 
-miR- 17 
-miR- 1451 
-miR- 2 Ob 
-miR-122-1 
-miR-128-2 
■miR- 9 -2 
-miR- 133b 
-miR-21 7 
-miR- 16c 
-miR-196-2 
miR-122-2 
-miR-15a 
-miR-21 8 -2 
-miR-1329 
miR-7-1 
-miR-218-1 



gcra-miR 
gga-miR 

gga-let 
gga-miR' 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-let 
gga-let 
gga-let 
gga-miR- 
gga-miR- 
gga-miR- 
gga-mlR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR 
gga-miR 
gga-miR 
gga-miR 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 
gga-miR- 



130e 
106 
27b 
7g 
;-7i 

146c 
:-16 
99a 
100 
181a 
23b 
16c 
454 
■1716 
184 
1677 
460 
181b 
199 
31 
-7c 
-7a 
7f 
10b 
218 
200a 
20b 
■10a 
200b 
205b 
203 
429 
1397 
33 
551 
217 
1329 
103 
■107 
■128 
-92 
1451 
24 
456 
1671 
383 
1698 
222 
489 
101 
130a 
1560 
18b 
148a 
15b 
15c 
1557 
1456 
18a 
1416 
30d 
30c 
133a 
122 
133b 
215 
206 
138 
193a 
la 
15a 
1684 
34b 



30b 
365 



Figure 4 miRNAs differentially expressed in the skeletal muscles of broilers and layers. A. Venn diagram demonstrating tine overlap of 
original sequenced reads between broilers and layers. B. Venn diagram demonstrating the overlap of sequenced reads with perfect genomic 
matches in broilers and layers. C Heat-map of miRNAs differentially expressed in broilers and layers based on the read counts obtained by deep 
sequencing. D. Heat-map of differentially expressed miRNAs confirmed by microarrays. Triplicate samples of total RNA from the skeletal muscle 
of ElO broilers and layers were used to perform miRNA microarray experiments. B, broiler chicken; L, layer chicken. 
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Figure 5 Validation of differentially expressed miRNAs using real-time RT-PCR. A. Real-time RT-PCR results for 17 miRNAs that were 
differentially expressed in broilers and layers. B. Relative abundance of 17 miRNAs in skeletal muscles of broilers and layers based on deep 
sequencing data. C Expression pattern of 17 miRNAs in skeletal muscles of broilers and layers based on microarray experiments. 
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Figure 6 Interaction network of differentially expressed miRNAs and their candidate targets. A. Workflow of interaction network analysis. 
Network construction can be divided into two components: miRNA-target interactions and target-target interactions. For miRNA-target 
interactions, candidate miRNA targets were predicted by TargetScan. Differentially expressed candidate targets were identified using mRNA 
microarrays that covered embryonic days 10, 12, 14 and 18. The final miRNA-target relations correspond to those mRNAs differentially expressed 
between broiler and layer that exhibited a pattern of expression opposite to that of the corresponding miRNAs. Target-target interaction pairs 
were extracted from the STRING database. A pairwise PCC was then calculated for each pair based on transcription profiles during skeletal 
muscle development to extract putative target-target interactions. B. The final integrated network. In the network, miRNAs are represented by 
yellow nodes and targets are represented by pink nodes. Blue lines denote miRNA-target interactions and red lines denote target-target 
interactions. 
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embryonic skeletal muscle that were differentially 
expressed between broilers and layers. In addition to 
analyzing mRNA transcriptomes on embryonic day 10 
(ElO), as was done for the miRNA transcriptome, 
embryos were analyzed on embryonic days 12, 14 and 
18 (E12, E14 and E18), yielding a total of 1057 non- 
redundant genes that were differentially expressed 
between broilers and layers (Additional file 12). To nar- 
row the field of candidate targets for the 24 miRNAs 
further, the analysis was restricted to those miRNA tar- 
gets that were differentially expressed between broilers 
and layers and had an expression pattern opposite to 
that of the corresponding miRNAs. Using these criteria, 
57 candidate targets for 16 miRNAs were identified 
(Additional file 13) and used for subsequent network 
analysis. 

In addition to interacting with miRNAs, targets should 
interact with each other. To investigate such interac- 
tions, protein-protein interactions (PPIs) of these puta- 
tive targets were extracted from the STRING database. 
Information concerning protein-protein interactions in 
bird species in the current PPIs database is very limited. 
Therefore, PPI data for human orthologs of these 
miRNA targets were utilized. It is widely accepted that 
some human PPIs may not be conserved in chickens, 
and protein interactions are time- and condition-speci- 
fic. To establish more reliable interactions, we applied a 
previously proposed referencing strategy in which PPIs 
are filtered with dynamic gene expression patterns 
[45,46]. The chicken gene expression dataset contained 
54 microarrays that covered nine developmental stages 
during skeletal muscle development of broiler and layer 
chickens (see Materials and Methods). To identify puta- 
tive target-target interactions, the Pearson Correlation 
Coefficient (PCC), which is known to provide informa- 
tion about the "shape" of gene expression changes 
[47,48], was used; an absolute PCC value of 0.3 was 
used as a cutoff. 

miRNA-target interactions and target-target interac- 
tions were integrated to construct possible regulatory 
networks (Figure 6B). One of the predicted miRNA-tar- 
get relationships presented in Figure 6B, between miR- 
27b and the target CYPIBI, has been reported pre- 
viously [49]. The remaining relationships are reported 
for the first time; therefore, this analysis predicts several 
candidates for future studies concerning miRNA-target 
function in controlling muscle development. In the pre- 
sented network, one major regulatory module involved 
13 miRNAs (yellow nodes) and 55 targets (pink nodes). 
Of these 13 miRNAs, five (miR-206, miR-la, miR499, 
miR-128 and miR-27b) have been reported to have a 
role during muscle development [30,50]. Little is known 
about the functional roles of the remaining eight (miR- 
31, miR-101, miR-200b, miR-lOb, miR-460, miR-15b, 



miR-16 and miR-203) during muscle development. 
However, analysis of their targets demonstrated that sev- 
eral were involved in myogenesis regulation, suggesting 
that these miRNAs could participate in regulating mus- 
cle development through their target genes. For exam- 
ple, miR-200b has eight predicted targets, among which 
are three genes RECK, SLC38A2 and Nr5a2, which 
encode proteins that are reported to be involved in mus- 
cle development [51-53]. 

It has been reported that activin A receptor type IIB 
(ACVR2B) plays an important role in regulating muscle 
development by interacting with a number of transform- 
ing growth factor-P (TGF-P) family members [54,55]. 
ACVR2B causes dramatic increases in muscle mass (up 
to 60% in two weeks) when injected into wild- type mice 
[29]. No miRNAs have been identified previously as reg- 
ulatory factors for ACVR2B, but the network analysis 
predicted that ACVR2B is a target of three miRNAs: 
gga-miR-101, gga-miR-la and gga-miR-499 (Figure 6B). 
It has been demonstrated that miR-1 is an important 
regulator of myogenesis [19,56]. miR-1 and ACVR2B 
had opposite expression patterns in skeletal muscle tis- 
sue from broiler and layer chickens (Figure 7B). There- 
fore, the target relationship between miR-1 and 
ACVR2B was validated using a luciferase reporter gene 
assay. As demonstrated in Figure 7C, the luciferase 
activity was significantly reduced when a miR-1 mimic 
was co-transfected with pGL3'ACVR2B'UTR containing 
a miR-1 targeting site into 293T cells, suggesting that 
miR-1 directly targets chicken ACVR2B UTR. Therefore, 
it is conceivable that miRNAs could be involved in regu- 
lating ACVR2B function in terms of controlling muscle 
development. Taken together, the results of the network 
analysis suggest that myogenesis is regulated by a com- 
plicated network, mediated by multiple miRNAs acting 
through the same target gene, and/or single miRNAs 
targeting multiple genes. 

Discussion 

Recent developments in high-throughput sequencing 
technology have enabled the miRNA transcriptome to 
be profiled in various organisms [36,57,58]. In the pre- 
sent study, Solexa deep sequencing was used to provide 
an extensive miRNA profile of the previously unexa- 
mined skeletal muscle of broiler and layer chicken lines. 
The sequence analysis identified 33 novel chicken miR- 
NAs and demonstrated that many miRNA precursors 
could generate multiple isoforms (isomiRs). Importantly, 
a comparison of miRNA transcriptomes allowed us to 
identify 16 known miRNAs and one novel miRNA that 
were differentially expressed in the skeletal muscles of 
broilers and layers. On the basis of the predicted targets 
of the 16 differentially expressed known miRNAs and 
eight muscle-related miRNAs, an interaction network 
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AGCTATGA < — mutant 
ACVR2B 3' UTR 5' ...CAAACUCAGUAUAUAACAUUCCA... 

I I I I I I I I 

gga-miR-1a 3' AUGUAUGAAGAAAUGUAAGGU 



ACVR2B 3' UTR 5' ...CAAACUCAGUAUAUAACAUUCCA.. 

I I I I I I I I 

gga-miR-1b 3' AUGUAUGAAGAAUUGUAAGGU 




miR-1 ACVR2B miR-l mimic - + - + - + 

Scramble + - + - + - 

Figure 7 miR-1 directly targets the chicl^en ACVR2B UTR. A. Schema of miR-1 binding site in cliicken AC\/R2B 3'-UTR sequence (seed 
sequence higliliglited in red). Mutated ACVR2B 3'-UTR eliminates tine seed binding site pointed by arrow. B. Expression of miR-1 and ACVR2B 
gene in sl<eletal muscle of broiler and layer chickens at embryonic day 18 were analyzed using real-time RT-PCR. C. Target validation using a 
luciferase reporter assay. 293T cells were co-transfected with miR-1 mimic or scramble double-strand small RNA and the reporter plasmid pGL3 
or pGL3-/\Cl//?28-UTR or pGL3-m/\Cl//?25-UTR, respectively. 



comprising these miRNAs and their candidate targets 
was constructed. 

The data presented in this report provide the first 
miRNA transcriptome profile of chicken skeletal muscle. 
Two hundred and thirty one known miRNAs and 29 
miRNA ''s were detected in skeletal muscles. gga-miR- 
206 was the most abundant miRNA in skeletal muscles 
of broilers (131,609 reads) and layers (222,998 reads), a 
result that is consistent with the well-established func- 
tion of miR-206 during skeletal muscle development 
[30]. Interestingly, Rathjen and colleagues recently per- 
formed miRNA profiling in chicken somites and demon- 
strated that among the 85 detectable known miRNAs, 
gga-miR-lOb was the most abundant (113,106 reads), 
whereas gga-miR-206 was much less abundant (259 
reads) [59]. Taken together, these observations suggest 
that miR-206 and miR-lOb could play important roles at 
different stages during muscle development. The 



expression level of myomiR miR-133 was lower than 
miR-206 in the skeletal muscle library, an outcome that 
could reflect differences in the roles of these miRNAs in 
terms of myogenesis regulation. Consistent with this 
interpretation, miR-206 has been shown to promote ske- 
letal muscle differentiation, whereas miR-133 regulates 
myogenesis by increasing muscle cell proliferation [19]. 

In addition to well-known myomiRs, recent studies 
have demonstrated that several other miRNAs are 
involved in regulating myogenesis. For example, a com- 
parison of miRNA expression profiles in proliferating 
myoblasts and differentiated myotubes revealed that 
miR-221 and miR-222 are down-regulated upon differ- 
entiation of primary and established myogenic cells, 
whereas miR-21, miR-103, miR-130, miR-99, miR-30 
and miR20 are up-regulated [19,33], suggesting that 
these miRNAs play important roles in the transition 
between proliferation and differentiation of muscle cells. 
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Interestingly, these same eight miRNAs were abundantly 
expressed in our sequencing libraries, indicating that 
they could play regulatory roles in controlling the differ- 
ence in skeletal muscle growth rates between broilers 
and layers during development. 

Seventeen miRNAs that were differentially expressed 
in the skeletal muscle of broiler and layer chickens were 
identified. Seven (miR-101, miR-lOa, miR-lOb, miR- 
1677, let-7f, miR-31, and miR-205b) were expressed at 
higher levels in layers, and ten (miR-203, miR-200b, 
miR-16c, miR-lSb, miR-15c, miR-460, miR-429, let-7c, 
miR-2188, and gga-miR-N2) were expressed at higher 
levels in broilers. Six of these miRNAs (miR-31, miR- 
10a, miR-lOb, miR-16C and two let-7 members) have 
been implicated in skeletal muscle regeneration or 
development [39-42]. Greco and colleagues demon- 
strated that miR-31 was induced in dystrophic (mdx) 
mice and in Duchenne muscular dystrophy patients, and 
in newborn mice and newly formed myofibers during 
postischemic regeneration, suggesting that it could be 
important in pathophysiological pathways that regulate 
muscle responses to damage and regeneration [42]. 
Recent studies have reported that miR-10 contributes to 
retinoic acid-induced smooth muscle cell differentiation 
[41], and may be important during the early stage of 
embryonic myogenesis [59]. Taken together, the 
approaches used in this study have identified a number 
of differentially expressed miRNAs that could exert 
novel functions in terms of regulating muscle cell prolif- 
eration and differentiation during development. Further 
investigations concerning the function of these miRNAs 
should facilitate our understanding of the regulatory 
roles of miRNAs in terms of controlling the divergent 
skeletal muscle growth rates of broiler and layer 
chickens. 

miRNAs exert their effects by interacting with target 
mRNAs. Therefore, target-predicting software (TargetS- 
can) was used to identify putative targets of these differ- 
entially expressed miRNAs and eight muscle-related 
miRNAs, before an interaction network of these miR- 
NAs and their candidate targets was constructed. The 
interaction networks predicted that ACVR2B is a target 
of gga-miR-101, gga-miR-la and gga-miR-499. Prior to 
this analysis, there have been no reports concerning 
associations between ACVR2B and miRNAs. The 
ACVR2B receptor signaling pathway mediates the func- 
tion of myostatin [60] and can regulate muscle growth 
in vivo [29]. ACVR2B haplotypes have been reported to 
be associated with muscle mass and strength in humans 
[60]. Furthermore, acute inhibition of myostatin/ 
ACVR2B signaling with the antagonist ACVR2B-Fc pre- 
serves skeletal muscle in mouse models of cancer 
cachexia [61]. ACVR2B was expressed at higher levels in 
the skeletal muscle of broilers than in layers at E18, 



indicating that ACVR2B could be related to the higher 
growth rate of broiler skeletal muscle. The results of 
this analysis indicate that the three putative miRNA reg- 
ulators of ACVR2B may be involved in this process. The 
analysis demonstrated that ACVR2B interacts with two 
other targets, CDR2 and GREMl. GERMl, a putative 
target of gga-miR-128, encodes a protein that is a BMP4 
antagonist and an effective regulator of myogenic pro- 
genitor proliferation [62]. 

RECK, the putative target of gga-miR-200b, is down- 
regulated by MyoD to facilitate myotube formation, and 
up-regulated by MRF4 to promote other aspects of 
myogenesis [51]. SLC38A2, encoding a sodium-coupled 
amino acid transporter, is the putative target of gga- 
miR-200b. SLC38A2 regulates proteolysis through phos- 
phoinositol 3-kinase, and provides a link among acido- 
sis, insulin resistance and protein wasting in skeletal 
muscle cells [52]. S0X8, the putative target of gga-miR- 
27b, acts as a specific negative regulator of skeletal mus- 
cle differentiation, possibly by interfering with the func- 
tion of myogenic basic helix-loop-helix proteins [63]. 
MEISl is the putative target gene of gga-miR-la and 
gga-miR-499. MEISl, together with PBXIA, facilitates 
binding of MyoD (a family of transcription factors with 
the remarkable ability to induce myogenesis in vitro and 
in vivo) to non-canonical E boxes in the myogenin gene 
to induce myogenesis [64]. Furthermore, a putative 
MEISl binding site is located in the minimal promoter 
of myostatin [65]. The CALDl gene, the putative target 
of gga-miR-27b, encodes two caldesmon-1 isoforms 
through alternative splicing: high molecular mass CaD 
(h-CaD), which is exclusively expressed in smooth mus- 
cle, and low molecular mass CaD (1-CaD), which is ubi- 
quitously expressed in all cell types except skeletal 
muscle. The h-CaD/l-CaD ratio can be used as a marker 
to monitor differentiating and pathological states of 
smooth muscles [66]. LATS2, a putative target of gga- 
miR-31, encodes a protein that has been reported to 
regulate the size of myocytes in the heart negatively 
[67]. NrSa2 is the putative target of gga-miR-200b, gga- 
miR-128 and gga-miR-27b. Its product, the nuclear 
receptor transcription factor Nr5a2, has been reported 
to function during skeletal muscle organization [53]. 
Therefore, although little is known about the specific 
functions of several of these miRNAs (e.g. miR-31, miR- 
101, miR-200b, miR-lOb, miR-460, miR-15b, miR-16 
and miR203) during muscle development, the close rela- 
tionship between their targets and myogenesis regula- 
tion demonstrates a potential role during muscle 
development. 

Given the significant functions of miRNAs in various 
biological processes, it is perhaps not surprising that 
miRNA biogenesis is tightly regulated at each stage of 
miRNA generation; in particular, a number of studies 
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have highlighted the complexity of post-transcriptional 
processing [9]. Sequence variations in mature miRNAs 
attributable to heterogeneity at the 5' and 3'-ends cre- 
ates an additional level of complexity in miRNA proces- 
sing [36]. A majority of miRNA genes have strand bias 
[68]. In some cases, miRNA genes have been found to 
generate similar amounts of miRNAs and their corre- 
sponding miRNA^'s [36]. The analysis of deep sequence 
tags identified several miRNAs that had read counts 
similar to those of their corresponding miRNA ''s, sug- 
gesting that these genes encode miRNAs on both arms 
of the precursor in skeletal muscle tissues. Although the 
functional significance of this observation has not been 
established experimentally, the fact that the miRNA and 
its miRNA'' are co-expressed at similar levels indicates 
that miRNAs serve independent functions in cells. 

This study identified 33 novel chicken miRNAs and 
analysis of the evolutionary conservation of these newly 
identified miRNAs revealed that only one is conserved 
in non- avian vertebrates and the remaining 32 are likely 
to be specific to bird and/or chicken lineages. Few 
newly identified miRNAs are conserved among verte- 
brates, whereas the majority of known miRNAs identi- 
fied using traditional cloning methods are abundantly 
expressed and relatively conserved during evolution 
[69-71]. Further support for this observation was pro- 
vided by a recent report concerning the identification of 
miRNAs in various organisms using high-throughput 
sequencing. This approach demonstrated that most 
newly identified miRNAs discovered using deep sequen- 
cing are present only in a small group of organisms 
[36]. Therefore, it is reasonable to hypothesize that 
these non-conserved miRNAs could play important 
roles in establishing and maintaining phenotypic diver- 
sity among different groups of organisms during evolu- 
tion. In addition, the bird and/or chicken lineage 
miRNAs reported in this study could have arisen during 
genetic selection, and function as key regulators of the 
differences in growth rates between broiler and layer 
chickens. A functional examination of novel and spe- 
cies-specific miRNAs is a challenge for future research 
and an important step in improving our understanding 
of the critical roles played by miRNAs during develop- 
ment and evolution. 

Conclusions 

The present study is the first to examine the chicken 
skeletal muscle miRNA transcriptome, and to evaluate 
miRNA function during skeletal muscle development 
through the identification of differentially expressed 
miRNAs between broiler and layer chickens, which have 
divergent skeletal muscle growth. Identification of novel 
miRNAs highlights the important function of low abun- 
dance and less conserved miRNAs during development 



of specific tissues. To investigate the functional roles of 
miRNAs during chicken skeletal muscle development, 
an interaction network of the differentially expressed 
miRNAs and their putative targets was constructed. 
This integrated analysis provides information that will 
aid further experimental investigations concerning miR- 
NAs and their targets during skeletal muscle 
development. 

Materials and methods 

Chicken embryo incubation and tissue collection 

Meat-type broiler eggs (Arbor Acres) and egg-type layer 
eggs (White Leghorn) were incubated at 37.5°C for 10 
or 18 days (ElO or E18). For Solexa sequencing and 
miRNA microarray analysis, skeletal muscles (pectoralis) 
were collected from broilers and layers at ElO; for 
mRNA microarray analysis, muscles were collected at 
ElO, E12, E14 and E18. Muscle samples were immedi- 
ately frozen in liquid nitrogen and stored at -80°C pend- 
ing RNA isolation. For analysis of the tissue expression 
pattern of novel miRNAs, tissues (brain, heart, liver, 
lung, breast muscle, intestine, kidney, fat, and stomach) 
were collected from layer chickens at E18. All embryo- 
nic manipulations were conducted in accordance with 
the protocols of the Chinese Academy of Medical 
Sciences and the Institutional Animal Care and Use 
Committee of Peking Union Medical College. 

Small RNA library construction and sequencing 

Total RNA was isolated from skeletal muscles using 
TRIzol reagent (Invitrogen) and precipitated overnight 
at -20°C. Approximately 20 (ig of total RNA from broi- 
ler and layer chickens was submitted to the Beijing 
Genomics Institute (BGI) for Solexa sequencing. In 
brief, sequencing was performed by fractionating total 
RNA using polyacrylamide gel electrophoresis (PAGE) 
to enrich for molecules in the range of 16-30 nt, and 
then ligated with proprietary adapters. Following adap- 
tor ligation, cDNA was synthesized from total RNA by 
reverse transcription and amplified with 15 PCR cycles 
to produce libraries for sequencing. 

Analysis of sequencing data 

After filtering low-quality reads and trimming the adap- 
tor sequences, totals of 2,700,003 and 2,576,562 reads 
were obtained for layers and broilers, respectively. 
Sequencing data were simplified by grouping identical 
sequence reads together, yielding 827,431 unique 
sequences. The unique sequence reads were mapped to 
the UCSC chicken genome galGal3 by ZOOM; 168,642 
uniq reads, corresponding to 3,541,220 sequences, were 
mapped to the genome with an exact match. 

The various types of ncRNAs or degradation products in 
the sequence library (Figure ID) were annotated by 
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reference to miRNAs from miRBase (version 16); coding 
exons based on RefSeq mRNA and repeat sequences based 
on RepeatMasker were obtained from UCSC; snoRNA, 
snRNA, rRNA and tRNA were based on Ensembl ncRNA 
data (version 54). Metazoan miRNA homologs of chicken 
miRNAs were identified in Anopheles gambiae, Ateles geof- 
froyi, Apis mellifera, Bombyx mori, Bos taurus, Caenorhab- 
ditis briggsae, Caenorhabditis elegans, Canis familiaris, 
Drosophila melanogaster, Danio rerio, Fugu rubripes, Gor- 
illa gorilla, Homo sapiens, Lemur catta, Lagothrix lagotri- 
cha, Monodelphis domestica, Macaca mulatta, Mus 
musculus, Macaca nemestrina, Ovis aries. Pan paniscus, 
Pongo pygmaeus. Pan troglodytes, Rattus norvegicus, Sagui- 
nus labiatus, Sus scrofa, Tetraodon nigroviridis, Xenopus 
laevis and Xenopus tropicalis. Sequencing reads represent- 
ing miRNA sequences often have untemplated nucleotides 
in the 3' end [72,73]. Therefore, miRNAs were annotated 
by identifying tags that were exactly matched to the 5' 19 
nt of known miRNAs. An analysis of the size distribution 
of these 5' end- matched sequences indicated that the most 
abundant size was 22 nt (Additional file 14). Those 20-25 
nt tags whose 5' 19 nt matched the 5' 19 nt of known miR- 
NAs were counted as copies of known miRNAs. The num- 
bers of chicken miRNA ''s, and miRNAs and miRNA ''s of 
other metazoans, were established using the same criterion. 
Other ncRNAs including snoRNA, snRNA, rRNA and 
tRNA, repeat sequences and mRNA degradation products 
were annotated on the basis of perfect matches. The anno- 
tation order was chicken miRNA, chicken miRNA'', 
metazoan miRNA homolog, metazoan miRNA'' homolog, 
snoRNA, snRNA, rRNA, tRNA, mRNA, and RepeatMas- 
ker. After each annotation step, only unmatched reads 
were used for the next annotation step. 

Simple perfect matches to miRBase (version 16) 
sequences were used to determine the 467 known 
chicken miRNAs and 77 miRNA^'s and their expression 
patterns (read numbers) in the data sets (Figure 2A and 
Additional file 5). For sequences that were not perfectly 
matched to known chicken miRNAs or miRNA^'s, 
metazoan miRNA homologs and metazoan miRNA'' 
homologs were identified using perfect sequence 
matches (Figure 2 A and Additional file 6). 

The deep sequencing data obtained were deposited in 
the GEO database with the accession number GSE20942. 

Prediction of novel miRNAs 

Novel miRNAs were identified using the miRDeep pack- 
age described by Friedlander et al., which can effectively 
distinguish miRNAs from other kinds of ncRNAs [37]. 
All sequences were mapped to the chicken genome 
(gal3) using megaBLAST, and only exactly matched 
sequences were retained for further analysis. As some 
miRNAs could lie in exonic regions of mRNAs [74,75], 
reads that aligned to more than five genomic positions 



or RepeatMasker annotation files were discarded. The 
remaining aligned reads were used as a reference, and 
potential precursor sequences were extracted from the 
chicken genome. The secondary structures of potential 
precursor sequences were predicted by RNAfold [76]. A 
FAST A file containing known mature metazoan miRNA 
sequences in miRBase was used as input to allow for 
conservation scoring. Using miRDeep, 222 putative miR- 
NAs were obtained, 189 of which mapped to known 
miRNAs in the chicken genome. The remaining 33 were 
novel chicken miRNAs (Table 3). Using the same cutoff 
on a permuted dataset, 16 putative miRNAs were 
obtained. Therefore, the corresponding signal-to-noise 
ratio was approximately 14:1 (222/16). 

Conservation analysis of miRNAs 

Genomic sequences for six vertebrates (hgl8, mmS, rn4, 
monDom4, xenTro2 and danRer4) were downloaded 
from the UCSC genome browser. BLASTN was used to 
identify regions of homology to chicken miRNA 
sequences in these genomes. From the BLAST results, 
sequences that covered more than 80% of the queried 
mature miRNA sequences and had fewer than two mis- 
matches in the covered region were selected. The seed 
regions of the miRNAs are more conserved, therefore 
covered regions were required to start from the 5' end. 
Finding "hit" sequences to mature miRNAs does not 
necessarily signify that the miRNAs are conserved as 
they may not be capable of forming hairpin structures. 
Accordingly, candidate sequences were extracted and all 
hairpin-like RNA structures encompassing small RNA 
sequence tags identified using RNAfold. Hairpin-like 
RNAs, whose mature miRNA regions were more than 
70% matched with the miRNA" regions, were accepted 
as homologs of chicken miRNAs. The conservation 
heat-map was constructed using Cluster 3.0 [77] and 
visualized in Tree View 1.60. In the heat-map (Additional 
file 2), the dark color represents 0, signifying that no 
homologs were identified in the corresponding species. 
The red color represents 1, which indicates the presence 
of homologs in the corresponding species. 

Identification of differentially expressed miRNAs on the 
basis of deep sequencing data 

miRNAs expressed at significantly different levels in 
broilers and layers were identified using the DEGseq 
package [38]. Using the likelihood ratio test model, pro- 
posed by Marioni et al. [78], and a cutoff of 1 x 10"^, 
102 known miRNAs were identified as being signifi- 
cantly differentially expressed. 

miRNA microarray 

Custom-designed miRNA microarrays were synthesized 
in situ by LC Sciences (Houston, USA) and used to 
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analyze miRNA expression patterns in the skeletal mus- 
cle of broilers and layers. Arrays contained 1721 DNA 
probes including a non-redundant set of probes comple- 
mentary to 440 known chicken miRNAs and 78 to 
known chicken miRNA'' s. Arrays included probes for 78 
predicted chicken snoRNAs, positive control probes for 
chicken U6 snRNA, 1124 unknown chicken small RNAs 
and negative controls for normalizing data with low- 
density signals. Hybridizations and scans were per- 
formed by LC Sciences; microarrays were scanned using 
an Axon GenePix 4000B Microarray Scanner. Data 
among arrays were normalized using a cyclic LOWESS 
(locally weighted regression) method [79] . The microar- 
ray data obtained were deposited in the GEO database 
with the accession number GSE20947. 

Quantitative real time RT-PCR 

Differentially expressed miRNAs identified using deep 
sequencing and microarrays were validated by stem-loop 
RT-PCR [80] using the stem-loop RT-PCR primers pre- 
sented in Additional file 15. Total RNA was isolated from 
skeletal muscles using TRIzol (Invitrogen), and genomic 
DNA contamination was removed by digesting with 
DNase I at 37°C for 30-40 min. DNase-treated RNAs 
were extracted using phenol/chloroform and precipitated 
with ethanol. Reverse transcriptase reactions contained 
RNA samples, 50 nM stem-loop RT primer, 1 x RT buf- 
fer, 0.25 mM each dNTPs, 3.33 U/ml MultiScribe reverse 
transcriptase and 0.25 U/ml RNase inhibitor. Reaction 
mixtures were incubated in a 9700 Thermocycler for 30 
min at 16°C, 30 min at 42°C, 5 min at 85°C, and held at 
4°C. Reverse transcriptase reactions including no-tem- 
plate controls and RT-minus controls were run in dupli- 
cate. Real-time PCR was performed using a standard 
SYBR Green PCR Master Mix (ABI). The reaction mix- 
tures were incubated in a 96-well plate at 95°C for 10 
min, followed by 40 cycles of 95°C for 15 s and 60°C for 
1 min. All reactions were run in triplicate. 

Network construction 

Network construction was divided into two components: 
miRNA-target interactions and target-target interactions. 
Starting with 16 known miRNAs validated as differen- 
tially expressed between broilers and layers, and eight 
muscle-related miRNAs (gga-miR-1, gga-miR-206, gga- 
miR-499, gga-miR-221, gga-miR-222, gga-miR-128, gga- 
miR-367 and gga-miR-27b), TargetScan (version 5.1) 
[43] was used to predict putative targets. Given the 
reported inverse correlation between miRNA and target 
expression patterns [44], the analysis was restricted to 
those differentially expressed miRNA targets whose 
mRNA expression pattern opposed that of the corre- 
sponding miRNAs. Expression levels of mRNA targets 
of miRNAs were measured using commercial Affymetrix 



Chicken Genome Arrays. Microarray experiments were 
carried out by CapitalBio Corporation (Beijing, China). 
Total RNA from skeletal muscles collected from broilers 
and layers at ElO, E12, E14 and E18, and prepared as 
described above, was hybridized in triplicate with three 
biological repeats from broilers and layers at each devel- 
opmental stage. Therefore, a total of 24 microarrays 
were used in the present study. Normalization was per- 
formed using RMA [81] software with a CDF file anno- 
tated by Dai et al. [82]. The Affymetrix GeneChip is a 
commonly used microarray platform for genome-wide 
expression studies. However, several genes/transcripts 
on the arrays are out of date owing to updates in gen- 
ome assemblies, causing problems when mapping the 
probes to new versions of the genome assembly [83]. To 
solve this problem, Dai et al. [82] aligned the probes to 
different sources of genome data to filter out proble- 
matic probes. The original. CEL files were re-annotated 
in the present study using annotations generated by Dai 
et al, ultimately obtaining 12,495 probe sets correspond- 
ing to 12,495 Entrez genes. Entrez genes whose expres- 
sion values changed more than 1.5 fold between broilers 
and layers were selected for each of four time points. To 
filter out genes that did not have identical expression 
profiles in each group, t-tests (p-value < 0.05) were used 
to obtain a final differentially expressed gene list. The 
microarray data obtained were deposited in the GEO 
database with the accession number GSE20990. 

In addition to interactions between miRNAs and their 
targets, interactions between miRNA targets were identi- 
fied. To establish target-target interactions, PPI data 
from the STRING database (version 8.0) were down- 
loaded. The chicken PPI data were limited, so human 
ortholog PPI data were used for these miRNA targets, 
applying a previously proposed referencing strategy that 
filtered the PPIs with gene expression patterns [45,46] 
to increase the reUability of the predicted interactions. 
The gene expression dataset used here contained nine 
time points: ElO, E12, E14, E18, Dayl (day of birth), W2 
(postnatal week 2), W4, W6 and W8 for broilers and 
layers. For each time point, there were three biological 
repeats for a total of 54 microarrays (Additional file 16). 
The generation of ElO, E12, E14 and E18 expression 
data is described above. Day 1, W2, W4, W6 and W8 
expression data were generated in a previous study [27]. 
For each target-target interaction pair extracted from 
the STRING database, a PCC value was calculated on 
the basis of the expression profiles from ElO to W8. An 
absolute PCC value of 0.3 (p-value < 0.05) was used as a 
cutoff to establish putative target-target interactions. 

Target validation using a luciferase reporter gene assay 

A pGL3-control vector (pGL3) was used for 3' UTR-luci- 
ferase reporter assays. The TargetScan Human database 
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http://www.targetscan.org/ was used to identify the pre- 
dicted miR-1 binding site. 3' UTR fragment of chicken 
ACVR2B containing a miR-1 binding site was ampUfied 
from chicken genomic DNA with primers [F: GCTCTA- 
GAGCTGGCCAGTTTTGAAGCAGAGGC (Xba I) and 
R: GCTCTAGAGCCCCCTGCTCACGGCTGTTGG (Xba 
I)] and cloned downstream of the luciferase gene to create 
the pGL3-luc-ACV7?2^-UTR constructs. miR-1 seed 
region mutations were generated by site-directed muta- 
genesis with primers (mut-F: CAAACTCAGTATA- 
TAAGCTATGAGTAAGGTTAGTATTGCAAAAC and 
mut-R: GCAATACTAACCTTACTCATAGCTTATA- 
TACTGAGTTTGATTGGT). Reporter assays were con- 
ducted in triphcate using 293T cells in 24-well plates. 
Transfections were performed with 150 ng of reporter 
plasmid and 50 ng of miR-1 mimic or scramble (Fugene; 
Roche). A pRL-TK reporter was used as an internal con- 
trol to normalize for transfection efficiencies. 

Additional material 
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